function VAR = doProxySVAR_biascorrected_single_EV(VAR)
% Assumes VAR and bias correction has already been estimated and is stored
% in structure=VAR.
% VAR.bet_bc = bias corrected betas
% Identify shocks using subsample with non-missing instrument

Xr      = VAR.X(VAR.I_m,:); % restricted subsample with non-missing instrument
Yr      = VAR.Y(VAR.I_m,:);
VAR.k   = 1;
VAR.mr  = VAR.m(VAR.I_m,:); 
[VAR.Tr,VAR.nr] = size(Yr);

% Residuals and Sigma for subsample with instruments 
% - uses bias corrected beta
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
VAR.resr   = Yr - [Xr ones(length(Xr),1)]*VAR.bet_bc;
VAR.Sigmar = (VAR.resr'*VAR.resr)/(VAR.Tr-VAR.nr*VAR.p-VAR.N);              % -VAR.N=30 for FE?

% Instrument
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Z = VAR.mr;

% Identification with ext instrument Z on RESTRICTED sample
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Phib    = [ones(length(Z),1) Z]\VAR.resr;           % Regress RF residuals on instrument and constant  
Phib    = Phib(2:end,:);                            % drop constant
VAR.Phib=Phib;
Phib11  = Phib(1:VAR.k,1:VAR.k);                    % First stage : HHD residual on instrument
Phib21  = Phib(1:VAR.k,VAR.k+1:VAR.Ny);             % Reduced form: lnY residual on instrument
b21ib11 = (Phib11\Phib21)';                         % RF/FS = s(21) / s(11)  
Sig11   = VAR.Sigmar(1:VAR.k,1:VAR.k);
Sig21   = VAR.Sigmar(VAR.k+1:VAR.Ny,1:VAR.k);
Sig22   = VAR.Sigmar(VAR.k+1:VAR.Ny,VAR.k+1:VAR.Ny);
ZZp     = b21ib11*Sig11*b21ib11'-(Sig21*b21ib11'+b21ib11*Sig21')+Sig22;    
b12b12p = (Sig21- b21ib11*Sig11)'*(ZZp\(Sig21- b21ib11*Sig11));            
b11b11p = Sig11-b12b12p;                                                  
b11     = sqrt(b11b11p);                                                 
VAR.b1  = [b11; b21ib11*b11];                                                                          
VAR.b21ib11=b21ib11;  


%Reliability
%  Sigmm = VAR.m'*VAR.m/VAR.T;
ED    = eye(VAR.k)*sum(sum(Z,2)~=0)/VAR.Tr;
mu1   = Z'*VAR.resr(:,1:VAR.k)/VAR.Tr;
 % PhiPhip = mu1*inv(b11b11p)*mu1';
 % VAR.RM    = inv(Sigmm)*PhiPhip*inv(ED)
Bi = [1/b11+(Sig21-Sig11*b21ib11)'*inv(ZZp)/b11*b21ib11 -(Sig21-Sig11*b21ib11)'*inv(ZZp)/b11];
VAR.et = (Bi*VAR.resr')';

PHI = mu1/b11;
GAM = inv(ED)*PHI;
E  = GAM*VAR.et(sum(Z,2)~=0);
V  = Z(sum(Z,2)~=0)-E;
VAR.RM = inv(E'*E+V'*V)*E'*E;


% Impulse Responses
%%%%%%%%%%%%%%%%%%%%
irs(VAR.p+1,:) = VAR.b1(:,1)/VAR.b1(1,1);
 for jj=2:VAR.irhor
 lvars = (irs(VAR.p+jj-1:-1:jj,:))';
 irs(VAR.p+jj,:) = lvars(:)'*VAR.bet_bc(1:VAR.p*VAR.Ny,:);     
 end
VAR.irs = irs(VAR.p+1:end,:); 


end



